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Abstract — Simulation of extracellular fields is one of the sub- 
stantial methods used in the area of computational neuroscience. 
Its most common usage is validation of experimental methods 
as EEG and extracellular spike recordings or modeling of 
physiological phenomena which can not be easily determined 
empirically. Continuous experimental work has been re-raising 
the importance of polarization effects between neuronal struc- 
tures to neuronal communication. As this effects rely on very 
small potential changes, better modeling methods are necessary 
to quantify the weak electrical fields in the microscopic scale in 
a more realistic way. 

An important factor of influence on local field effects in the 
hippocampal formation is the heterogeneous resistivity of extra- 
cellular tissue. The vast majority of modeling studies consider the 
extracellular space to be homogeneous while experimentally, it 
has been shown that the stratum pyramidale has two times higher 
resistivity than other hippocampal layers. Common simulation 
methods for extracellular electrical fields based on the point 
source approximation are bound to describe the resistance of 
the space with a single scalar. We propose that models should 
be based on the space- and time-dependent Maxwell equations 
(Partial Differential Equations, PDEs) in order to account for 
heterogeneous properties of the extracellular space and specific 
arrangements of neurons in dense hippocampal layers. 

To demonstrate the influence of heterogeneous extracellular 
resistivity and neuronal spatial orientation on modeling results, 
we combine solutions of classical compartment models with 
spatiotemporal PDEs solved by the Finite Element Method 
(FEM). With the help of these methods, we show that the 
inclusion of heterogeneous resistivity has a substantial impact 
on voltages in close proximity to emitting hippocampal neurons, 
substantially increasing the change in extracellular potentials 
compared to the homogeneous variant. 

Index Terms — extracellular fields finite element method neu- 
ronal arrangement 



I. Introduction 

Numerous computational studies have investigated time- 
varying currents in homogeneous extracellular space (D, lEI, 
(J I as well as the role of neuronal morphology in uniform 
electric field stimulation O, O, 0. However, most studies 
use models of the extracellular milieu that may not be ac- 
curately applied when modeling brain structures with highly 
diverse extracellular resistivity and neuronal arrangement. 



As emphasized by Lopez- Aguado and Bokil f?!, fS^ it 
has been traditionally neglected that currents propagate in all 
directions in an extracellular medium and that inward and 
outward currents originate from tissue regions having large re- 
sistivity differences. A usual argument for this approach is that 
resistivity influences extracellular electrical fields minimally 
and that extracellular space can be assumed as homogeneous 
in current- source density (CSD). However, extracellular non- 
homogeneous resistivity has been shown experimentally in 
several regions of the brain, for example, in the hippocampus 
and the cerebellum |9|. Additionally, tissue swelling has 
been observed after intense neural activity that in turn could 
lead to an increase in extracellular resistivity [10|. Previous 
studies [T|, ifTTIl . (T2I, (13 1 have attempted to examine how 
extracellular electrical fields affects neuronal activity although 
with the help of quasi-static approximation and an assumed 
homogeneous tissue resistivity. 

To quantify the effect of extracellular heterogeneous resis- 
tivity and neuronal spatial distribution on strength of neu- 
ral fields, we are proposing a simple modeling pathway to 
couple compartment-based neural models with the COMSOL 
Multiphysics simulation environment. Here, we solve time- 
dependent Maxwell's equations using the FEM to analyze the 
change of electrical fields as it occurs in the extracellular space 
surrounding neurons. 

Our results indicate that inhomogeneous resistivity of the 
extracellular milieu significantly influences the change of ex- 
tracellular potentials (EPs) in the hippocampus. By computing 
the resulting voltage change due to outgoing transmembrane 
currents of an exemplar CAl pyramidal cell model in ho- 
mogeneous and inhomogeneous extracellular resistivity, we 
observed a maximal difference in EP change of 60% in the 
hippocampal pyramidal layer. Furthermore, the here proposed 
method offers the possibility to efficiently simulate the effects 
of superimposed extracellular potentials created by neurons in 
divergent positions relative to each other. We will also discuss 
advantages and drawbacks of the used method and propose 
alternatives and possible improvements towards more realistic 
modeling of electrical fields of the brain. 
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II. Materials and Methods 

A. Pyramidal neuron model 

To demonstrate the effect of trans-membrane currents on 
the effect of extracellular fields, a Hodgkin-Huxley (HH)-like 
hippocampal CAl pyramidal cells was adapted from [14] with 
the addition of a current Ih from ifTSl and a current Im from 

m. 
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Ih current: 

iXoo = [l+exp((F + 76)/7)]-\ 
Tu = 10V[237exp((V + 50)/12)+ 

17exp(-(V + 50)/25)] + 0.6, 
h{t,V) = ghu{t,V) {V-Eh). 

Im current: 

= [1 + exp(-(y + 22.53)/10)]-\ 
Ts = 4135.7/[164.64exp((V - 0.05) • 0.12)+ 

0.33exp(-(V - 0.05)/10)] + 35.66, 

Im{t,V) = gms{t,V) (V-Ek). 
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Our model contained one somatic (r = 10 /im), 20 dendritic 
(r = 1-3 /im, 1 = 5 /im), two axon initial segments, AISi 
(r = 0.5 /im, 1 = 60 jim ), AIS2 (r = 0.5 /im, 1 = 60 /im), 
and finally 32 axonal compartments (r = 0.5 /im, 1 = 5 /im), 
see Figure [T] Note that here AISi refers to the part of AIS 
from to 60 /im, and AIS2 from 60 to 120 /im from the cell 
body. Additionally, AIS2 compartments contain high Na'^ 
channel density (O- 20 dendritic compartments are used for 
the dendritic tree, while additional compartments are included 
in the branching analysis. The Na^ channel density was 
varied between 986 and 2943 (corresponding to AIS2). 
The following conductance values were used: gNa = 8, gkdr 
= = I mS/cm?' for the soma; gca = 10, gxc - 

QKaHP = 0.8, go = 0.625, gNa = 0.07, gn = 0.2 mS/cm^ 
for the dendrite; gNd = 50, gkdr = 10, gm = 10 mS/cm?' for 
the AIS2 and gNa = 9, gkdr = 10, = 10 mS/cm?' for the 
rest of the axon compartments. The leakage conductance was 
set to 0.1 for all of the compartments. Finally, the equilibrium 
potentials were set to ENa = 60 mV, Ek = -85 mV, E^ = -43 
mV and Eieakage = -65 mV. 




Fig. 1. Mapping trans-membrane currents from compartmental model to 
FEM boundary mesh. A compartmental model containing 20 dendritic, one 
soma, two axon initial segments and 32 axonal compartments. Ionic currents 
for each compartment are shown on the left. The compartmental model outputs 
I{t) which is used as current boundary sources Qj{t) in the corresponding 
three-dimensional volumes in the FEM model. 



B. Creation of morphology and model coupling 

The three-dimensional neuronal geometry was constructed 
in COMSOL Multiphysics 4.3 with the help of the interface 
to MATLAB ("LiveLink") by morphological additions and 
boolean unification of simple geometric volumes. As we aim 
to represent the 3D-morphology as an exact counterpart of the 
compartmental model, each section is recreated as a cylinder 
with the same length and diameter as in the compartmental 
definition. If the cylinders were added on top of each other 
with no change of the rotation vector, as shown in Figure [T] no 
joining geometrical primitives were added in between them. 
If the rotation differs, as for example in Figure |4] a sphere 
is added in between the cylinders, followed by a removal of 
the interior boundaries of both the cylinders and the sphere. 
The reason for constructing the geometry in this way is that 
the mesh engine otherwise respects the internal boundaries 
such that the resulting mesh becomes unnecessary complex. 
In some cases the occurrence of multiple internal boundaries 
can even make the meshing procedure abort. 

Ionic currents of single neuronal compartments model were 
determined by solving Hodgkin-Huxley Ordinary Differential 
Equations (HH ODEs) specified in equations 1 to 6 using a 



Runge-Kutta algorithm from the MATLAB ODE-suite (Math- 
works) with a constant step- size T. The sum of all currents 
per compartment for each time step of the simulation /q is 
afterwards normalized to an absolute value of units A/m? 
and stored in a matrix / of size N compartments x Numesteps- 

Next, by using the interface to COMSOL, each row of 
the matrix / was mapped adequately to the corresponding 
cylindrical domain as a boundary current source Qj{t). Note 
that we hereby assume that transmembrane currents are the 
only cause of change of extracellular potential, which is surely 
not the case in a real neuron, as for example synaptic calcium- 
mediated currents are suspected to contribute to a large fraction 
of the extracellular signature iflSl . 

The extracellular volume was modeled by cylindrical ob- 
jects covering the neuronal morphology with either homo- 
geneous (0.3 S/m) or heterogeneous resistivity as shown in 
Figure [2]\. In this case we constructed the extracellular volume 
by taking the union of cylindrical objects with increasing 
resistivity in the ?/-axis, according to |7|. It is assumed that the 
conductivity of extracellular tissue is frequency-independent in 
the used range of neural activity (10-100 Hz) [|T9l . 
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C. Electrostatic formulation and the Finite Element Method 

To simulate non-homogeneous distributions of electrical 
fields produced by single neurons and neuronal networks, 
we used an electrostatic formulation of Maxwell's equations 
discretized with finite elements. Chiefly, in the finite element 
approach. Maxwell's equations are solved by discretizing the 
incorporated volumes (in this case the neuronal compartments) 
into finite tetrahedral volume elements [20|. We seek the elec- 
tric field intensity E in terms of the electric scalar potential 

E = -VK (7) 

The relevant dynamic form of the continuity equation with 
current sources is given by 

V-J = -|^+Q„ (8) 

with J and p the current density and electric charge density, 
respectively. Further constitutive relations include 

D = eo^rE, (9) 

and Ohm's law 

J = crE, (10) 

in which D denotes the electric flux density. Finally, Gauss' 
law states that 

V-D = p. (11) 

Upon taking the divergence of JTOl) and using the continuity 
equation (O we get 

V-aE = -|^ + Q,. (12) 

Rewriting the electric charge density using Gauss' law together 
with the constitutive relation (O and finally applying the gauge 
condition © twice we arrive at the time-dependent potential 
formulation 

- V • (^a\/V + eoSr^WV^ = Qj. (13) 

This is the formulation used in COMSOL Multiphysics (211. 
The values for the electric conductance a and the relative 
permittivity Sr were obtained from [2] . The source currents Qj 
were computed from the compartmental model as described 
above. 

The formulation ([T3]) is efficiently solved by COMSOL' s 
"Time discrete solver", which is based on the observation 
that the variable W := AV satisfies a simple ODE. Solving 
for W in an independent manner up to time t, it is then 
straightforward to solve a single static PDF to arrive at the 
potential V itself. 

As for boundary conditions we took homogeneous Neu- 
mann conditions (electric isolation) everywhere except for in 
a single point which we choosed to be ground (V = 0). 
In all our simulations this point was placed at the axis of 
rotation of the enclosing cylindrical extracellular space, and 
underneath the neuronal geometry. This procedure ensures that 
the formulation has a unique solution (it is otherwise only 
specified up to a constant). 



A tetrahedral mesh was applied to discretize space (using 
the "finest" mesh setting; resolution of curvature 0.2, reso- 
lution of narrow regions 0.8). The simulations were verified 
against coarser mesh settings in order to ensure a practically 
converged solution. As a final note, in the Time discrete solver 
the time step was set to the same step size T as used in 
the ODF-based solution of the Hodgkin-Huxley equations, 
thus ensuring a correct transition between both simulation 
environments. 



III. Results 

A. The effect of heterogeneous extracellular fields 

In order to examine the influence of heterogeneous ex- 
tracellular space, we first constructed a three-dimensional 
active neuron model in an homogeneous and heterogeneous 
extracellular milieu (Figure O. Neurons in the hippocampus 
have an intricate spatial orientation that propitiates strong field 
potentials: high density of pyramidal cell dendrites running 
in parallel in the stratum radiatum (SR), densely packed 
pyramidal cell somas in the stratum pyramidale (SP) while 
pyramidal cell axons run almost in parallel or crossing each 
other in SP or stratum oriens (SO). 

We measured the voltage on four defined point probes 
placed parallel to the dendrites, soma, AIS2 and axon terminal 
compartments during the peak of an action potential (AP) of 
80m V, by varying the distance between active neuron and point 
probes from 1 pm to 80 pm (Figure [JK). 

In the first set of simulations, we analyzed the effect of 
the aforementioned four neuronal regions on the defined point 
probes assuming a widely accepted homogeneous resistivity of 
350 Q cm [2J, |3|. Note that here applied boundary currents 
of three-dimensional neuronal compartments correspond to the 
peak transmembrane current during an AP. By doing so, the 
peak voltage of 0.25 mV was obtained in the point probe 
parallel to the AIS2 compartment, followed by soma, dendrite 
and axon terminal. Consequently, the heterogeneous case was 
examined by placing the neuron and point probes into a 
heterogeneous extracellular space representing hippocampal 
spatial order, where resistivity values for different strata were 
obtained from |7| and shown in Figure |2]\. The active neuron 
was positioned in the center of SP and point probes were 
moved along the x axis of the extracellular space. In the case 
of non-homogeneous extracellular resistivity (Figure |2]l!) the 
largest voltage change was measured in point probes parallel 
to AIS2, although the values in point probes placed in SP 
were 60% larger in close distance than in the homogeneous 
extracellular space scenario and 28% higher considering the 
spatial mean over the total distance. The point probes placed in 
SO and SR were as well affected by the higher restivity of the 
pyramidal layer, showing an average increase of 4% in parallel 
to the axon and 7% in parallel to the dendrite. The appreciable 
difference between the voltage changes suggests that non- 
homogeneous distribution of resistivity is an important aspect 
of extracellular field effects. 
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Fig. 2. Heterogeneous extracellular resistivity modulates the strength 
of polarization effects in proximity of the firing neuron. A, Schematic 
representation of a model neuron impact on defined point probes. Point probes 
act as receivers and are being shifted away from the active neuron along 
the X-axis (left). Non-homogeneous resistivity distribution of hippocampal 
regions (right). B, Relationship between membrane potential at different point 
probes affected by dendritic (D), somatic (S), axon initial segment (A/S'2) and 
axon terminal (A) neuronal compartments (color coded as in A) for distance 
varying between and 80 /im. Dashed: the resulting voltage change assuming 
homogenous extracellular resistivity (= 350 fl cm), while solid represents the 
result of assuming a heterogenous resistivity. 



B. Analysis in respect to spatial orientation and neural mor- 
phology 

According to the superposition principle electrical fields 
(EFs) are produced by summation of single neuronal activity. 
Thus, the mutual interaction between neuron and field is 
strongly modulated by the spatial orientation of neuronal 
assemblies ll22l . To analyze EFs considering a realistic spatial 
order, we began simulating two neurons with parallel axons 
(Figure [3]\) where NeuroUa represents an active structure 
with boundary current source Qj mapped to its surface while 
Neuroup is a passive measurement structure. We measured 
Ve at four points along the axon of NeuroUp and detected 
the maximal voltage amplitude of 0.28 mV. Note that here 
Ve corresponds to extracellular voltage measured on the cell 
membrane. Next we tested the influence of four neighboring 
cells with parallel axons on the NeuroUp (Figure [3^). In this 
case, having four active neuronal neighbors (distance between 
the axons =10 jim) firing non- synchronously, a maximal 
voltage of 0.49 mV in homogeneous and 0.68 mV in heteroge- 
neous extracellular medium of NeuroUp was computed. Syn- 
chronous firing of adjacent neurons produced a peak voltage 
of 0.62 mV in homogeneous and 0.84 mV in heterogeneous 
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Fig. 3. Analysis of extracellular fields in dependence of spike timing 
and spatial orientation. A. Potential change calculated at different regions 
(colored circles) of a passive neuron (Neuroup) in response to EPs triggered 
from action potentials in the initial phase of an active neuron. Lines in the 
right panel are color coded as the circles in the left panel. Somas were 
placed at the stratum pyramidale and extracellular resistivity followed the 
same distribution as shown in Figure 2C. Upper trace of the right panel 
shows potential changes in homogeneous, while lower trace in heterogeneous 
extracellular space. B, Same as in A but the Neuroup surrounded by four 
active neurons. Potential change at different regions of Neuroup when four 
active neurons fire asynchronously (middle panel) and synchronously (right 
panel). Insets on the middle and right panels show the firing of active neurons). 
C, Same as in B, but the active neurons where oblique (but not intersecting, 
nearest interaxonal distance = 2 iim) at the position of AIS2 to Neuroup. 



medium of AIS2 • When axons crossed each other at the AIS2 
region (Figure [3]^), we observed a maximal voltage amplitude 
of 0.89 mV in the AIS2 of Axorip during asynchronous 
activity of the neighboring cells, whereas synchronous activity 
produced 1.41 mV in heterogeneous medium. These results 
suggest that for analyzed neuronal arrangement soma and 
AIS2 have considerable influence on the strength of neural 
fields for small distances (< 10 /im) and that values in 
heterogeneous extracellular space are in mean ~ 25% larger 
than in homogeneous extracellular medium when averaged 
over time. 

Furthermore, we simulated the influence of pyramidal 
neuron position within different hippocampal layers on the 
strength of neural fields (Figure |4]l. We positioned somas in 
SP with the smallest distance between the two axon initial 
segments confined to SO. This geometrical arrangement was 
motivated by confocal images of two proximal pyramidal 
neurons located parallel to each other with an axonal bend 
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Fig. 4. Modeling the effect of bended axons. A, Confocal image of 
two neurons in the hippocampus filled with neurobiotin for post hoc spatial 
analysis showing two parallel neurons with bent axons situated at a minimal 
distance of 2 iim. B, A simplified representation of both neurons when the 
nearest interaxonal region was confined in stratum oriens. In this case, an 
action potential in the active neuron generates a voltage transient of 0.25 
mV at the passive neuron soma and 0.35 mV at the neuron AIS2. The figure 
shows the extracellular voltage distribution at the peak amplitude of the active 
cell action potential. C, Same as B, but axon initial segments were confined at 
the stratum pyramidale (higher resistivity). In this case, a 0.25 mV and 0.45 
mV voltage transients were observed at the passive neuron soma and axon, 
respectively. 

Starting about 50 /im behind the soma (Figure |4]\). As 
known from previous analytical studies [23 1 the bend of 
axonal structures generally amplifies endogenous fields around 
neurons. We were able to confirm this effect in our simulation 
as we simulated the field effect resulting from trans-membrane 
currents of an active neuron measured on the parallel neuron. 
In this scenario, peak voltage of 0.35 mV was calculated at the 
position of AIS2 while the smallest voltages were registered 
at the passive neurons between the dendrites (0.12 mV) in SR 
and axon terminals (0.16 mV). If AIS2 were confined to SP, 
the potential computed at the passive neuron was 25% larger 
than in the case in which AIS2 were positioned in the SO 
(Figured^). 

IV. Discussion 

In this work, we use the Finite Element Method coupled 
to the HH-equations to simulate how neuronal geometry, 
arrangement and heterogeneous extracellular properties affect 
the strengths of neural fields. We first show that there is a 
difference in voltage transients produced by firing of neigh- 
boring cells if the extracellular space is considered to be 
homogeneous or heterogeneous. Additionally, we demonstrate 
that the spatial orientation of specific cellular compartments is 
an important determinant of the strength of neural fields. 



In our computations, the highest change of extracellular 
potential arises in the pyramidal layer, in proximity to the 
AIS2 compartment. Action potentials are generally initiated in 
the AIS due to a higher density of voltage-gated Na~^ channels 
C3, II24II . |l25l, which is reflected by the parameter settings 
for this compartment in the occupied neuronal model. Ad- 
ditionally, the spatial arrangement of hippocampal pyramidal 
neurons also propitiates the proximity of AIS2 in both SP and 
SO (see Figure IS |26|). 

Holt and Koch 1 27 ] showed that interactions near cell bodies 
are more important than interactions between axons by using 
standard one-dimensional cable theory and volume conductor 
theory. Another study [|28l reported a 4.5 mV change (in the 
AIS) caused by extracellular interactions, a change more than 
4x larger than what we have found in our simulations of 
heterogeneous extracellular media. Additionally, by applying 
analytical methods, Bokil et al |8| have shown that, in the 
olfactory system, an action potential of 100 mV amplitude 
in one axon could produce depolarization in other axons in 
the bundle sufficient to initiate an action potential. However, 
using FEM simulations it was not possible to trigger spikes in 
neurons solely by electrical fields mediated by in- and outflow 
of trans-membrane currents, in agreement with the work by 
Traub and colleagues (in which maximal voltages in a sink 
axon during synchronized activity of four neighboring neurons 
is ^ 1.2mV ). 

Hence, simulations relying on the point source approxima- 
tion to describe hippocampal neural fields may be distorted as 
the potential change is more than 28% greater in stratum pyra- 
midale and in average 6% greater in other hippocampal layers 
if a heterogeneous tissue is used instead of a homogeneous 

one m, cn. 

For varying extracellular resistivity a numerical procedure 
referred to as 'the method of images' [|29l has been proposed 
as an extension to the point source equation. Although this 
method has been used in computations of extracellular action 
potentials in a previous study [HJ, its practical use is limited 
to a rather low number of resistivity layers and non-complex 
geometry. 

The requirement of FEM to model passive current flow 
between neurons was suggested elsewhere fS], f3l, but im- 
plementation issues and the lack of adequate software tools 
may have precluded its usage in the past. 

Extracellular resistivity could also contribute to the ampli- 
tude of local field potentials (LFP), and in fact, in vitro LFP 
registered in hippocampal slices are greater in the SP than in 
other layers [30|. Interestingly, we observed that hippocampal 
slices in interface-type recording chambers (where slices are 
not completely submerged in the buffer solution [1311 . [|32l ) are 
more than 60% less conductive than in the chambers where 
slices are submerged (unpublished results). This could help 
explaining why LFP recorded in interface-type cambers are 
far greater than in submerged-type counter parts [16], [32] . 

Here we show that the use of FEM software (COMSOL 
Multyphysics), with an interface to Hodgkin and Huxley- 
type ODE model in Matlab, is a powerful tool to verify 
macroscopic effects of neural fields. However, this approach 
may fail to simulate small nuances of extracellular interactions 
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(e.g. ion channels apposing active compartments may suffer 
more from the effect of passive current flow than channels 
in the opposite side). Nonetheless, simulations using PDEs 
could increase the level of realistic models of membrane 
dynamics. The idea of translating membrane dynamics to 
PDE was proposed by Hodgkin and Huxley themselves |[33l ^ 
and later pursued by Kashef and Bellman |34|. However, 
implementing membrane dynamics in combination with the 
Maxwell equation interface of the PDE simulator has been, to 
our experience, quite cumbersome. For example, due to the 
relatively high computational cost of the solution phase of the 
finite element method, a representation of fully reconstructed 
morphologies with volumetric cylindrical elements was not 
possible. The high amount of (small) cylindrical elements 
required to accurately model complex dendritic branches typ- 
ically caused the meshing algorithm to break. 

The software used in our simulations (COMSOL Multi- 
physics) has helped to popularize the use of FEM in neurosci- 
entific problems ll35l . ll36ll . However, the software is geared 
towards industrial applications and the steep learning curve 
associated to the adaptation of COMSOL to neuroscience 
related problems may preclude its widespread usage by the 
neuroscience community. Possible combination of Neuron 
Simulation Environment (NSE) with a FEM simulator would 
offer great possibilities by allowing both realistic neural mor- 
phologies and extracellular space modeling. However, cur- 
rently the existence of a direct interface between NSE and a 
FEM simulator has not been reported; future studies should 
focus on this task. Furthermore, hippocampal extracellular 
stimulation has been proposed as therapeutic strategy to con- 
trol several disorders, including temporal lobe epilepsy [l37l . 
lISSl . Developing models that consider as many as possible 
physical tissue properties would help to improve extracellular 
stimulation in epilepsy patients. In summary, our work adds to 
the recently published studies attempting to reveal important 
parameters determining the strength of extracellular electric 
fields. Models that do not use space and time-dependent 
differential equations when modeling neuronal interactions 
may have failed to replicate the changes in measured voltages 
caused by passive current flow in heterogeneous extracellular 
tissue, especially when more than two neurons are modeled 
simultaneously. 
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